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Abstract One main issue, when numerically integrating autonomous Hamiltonian 
systems, is the long-term conservation of some of its invariants, among which the 
Hamiltonian function itself. Recently, a new class of methods, named Hamiltonian 
Boundary Value Methods (HBVMs) has been introduced and analysed [5|, which are 
able to exactly preserve polynomial Hamiltonians of arbitrarily high degree. We here 
study a further property of such methods, namely that of having, when cast as Runge- 
Kutta methods, a matrix of the Butcher tableau with the same spectrum (apart the zero 
eigenvalues) as that of the corresponding Gauss-Legendre method, independently of 
the considered abscissae. Consequently, HBVMs are always perfectly A-stable meth- 
ods. Moreover, this allows their efficient blended implementation, for solving the 
generated discrete problems. 
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1 Introduction 

Hamiltonian problems are of great interest in many fields of application, ranging from 
the macro-scale of celestial mechanics, to the micro-scale of molecular dynamics. 
They have been deeply studied, from the point of view of the mathematical analysis, 
since two centuries. Their numerical solution is a more recent field of investigation, 
which has led to define symplectic methods, i.e., the simplecticity of the discrete 
map, considering that, for the continuous flow, simplecticity implies the conservation 
of H(y). However, the conservation of the Hamiltonian and simplecticity of the flow 
cannot be satisfied at the same time unless the integrator produces the exact solution 
(see lfl6l page 379]). More recently, the conservation of energy has been approached 
by means of the definition of the discrete line integral, in a series of papers 118111911201 
|2T1|221 . leading to the definition of Hamiltonian Boundary Value Methods (HBVMs) 
[3]|4][5]|6), which are a class of methods able to preserve, for the discrete solution, 
polynomial Hamiltonians of arbitrarily high degree (and then, a practical conserva- 
tion of any sufficiently differentiable Hamiltonian). In more details, in Q, HBVMs 
based on Lobatto nodes have been analysed, whereas in [6] HBVMs based on Gauss- 
Legendre abscissae have been considered. In the last reference, it has been actually 
shown that both formulae are essentially equivalent to each other, in the sense that 
they share the same order (twice the number of fundamental stages) and stability 
properties, and both methods provide the very same numerical solution, when the 
number of the so called silent stages tends to infitinyQ In this paper this conclusion 
if further supported, since we prove that all such methods, when cast as Runge-Kutta 
methods, have the corresponding matrix of the tableau, whose nonzero eigenvalues 
coincide with those of the corresponding Gauss-Legendre formula (isospectral prop- 
erty of HBVMs). 

This property can be used to define an efficient iteration for solving the discrete 
problems generated by the methods, via their blended implementation. Indeed, after 
posing HBVMs in block BVM form, they can be recast in the framework of blended 
implicit methods, which have been studied in a series of papers |[23l7l8ll9l [TOll 1 1 II 1 2l 
fT3in~5l (see also C. Magherini's PhD Thesis [23]). The latter methods have been suc- 
cessfully implemented in the two computational codes BiM and BiMD [241; the latter 
code is also included in the current release of the "Test Set for I VP Solvers" l25l . 

With this premise, the structure of the paper is the following: in Section[2]the basic 
facts about HBVMs are recalled; in Section [3] we state the main result of this paper, 
concerning the isospectral property; in Section|4]the discrete problem to be actually 
solved is defined; in Section|5]it is shown that a corresponding blended iteration can 
be devised for its efficient solution, which can be tuned by choosing a free parameter; 
in Section [6] the optimal choice of the free parameter is done, on the basis of the 
isospectral property of HBVM(A:, s), by using a linear analysis of convergence; finally, 
in Section|7]a few concluding remarks are given. 



Actually, they both provide the same numerical solution also when the Hamiltonian is a polynomial 
and the number of silent stages is high enough to ensure the conservation property of the Hamiltonian 
function itself (see )2.81 ). 
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2 Hamiltonian Boundary Value Methods 

The arguments in this section are worked out starting from the arguments used in ||5] 
|6l to introduce and analyse HBVMs. We consider canonical Hamiltonian problems 
in the form 

y = JVH(y), y(t )=yoem. 2m , (2.1) 

where / is a skew-symmetric constant matrix, and the Hamiltonian H(y) is assumed 
to be sufficiently differentiable. The key formula which HBVMs rely on, is the line 
integral and the related property of conservative vector fields: 

H(y 1 )-H(y )=h C 6{t + zh) T VH(a{t + Th))dz, (2.2) 
Jo 

for any y\ £ M 2 '", where a is any smooth function such that 

ff(/o)=yo, a{t +h)=yi. (2.3) 

Here we consider the case where a(t) is a polynomial of degree s, yielding an ap- 
proximation to the true solution y(f) in the time interval [to, to +h]. The numerical 
approximation for the subsequent time-step, yi, is then defined by ( 12.3b . 
After introducing a set of s distinct abscissae 

0<ci,...,c,<l, (2.4) 

we set 

Yi = a(t + Cih), i=l,...,s, (2.5) 

so that o(t) may be thought of as an interpolation polynomial, interpolating the fun- 
damental stages Yi, i= \,...,s, at the abscissae ( 12.41 i. We observe that, due to ( 12.31 ), 
a(t) also interpolates the initial condition yo. 

Remark 1 Sometimes, the interpolation at to is explicitly required. In such a case, 
the extra abscissa cq — is formally added to \2.4\ . This is the case, for example, of 
a Lobatto distribution of the abscissae ^5]/. 

Let us consider the following expansions of o(t) and a(t) for t € [to,to + h]: 

S s „ T 

CT(fo + TA) = £y;P,-(T), a{to + xh)=yo + hY,7j Pj{x)6x, (2.6) 

where {Pj(t)} is a suitable basis of the vector space of polynomials of degree at most 
s — 1 and the (vector) coefficients {jj} are to be determined. We shall consider an 
orthonormal basis of polynomials on the interval [0, 1], i.e.: 

[ Pi(t)Pj(t)dt = 8 u , i.j 1 v. (2.7) 

Jo 

where 5 !; is the Kronecker symbol, and P, (f) has degree i— 1. Such a basis can be 
readily obtained as 

F,(t) = y/2i-lPi-i(t), i = 1 J, 

with Pi-i(t) the shifted Legendre polynomial, of degree i — 1, on the interval [0, 1]. 
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Remark 2 From the properties of shifted Legendre polynomials (see, e.g., fiTj or the 
Appendix in ^J), one readily obtains that the polynomials {Pj(t)} satisfy the three- 
terms recurrence relation: 

Pi(t) = l, ft(/) = >/3(2f-l), 

We shall also assume that H(y) is a polynomial, which implies that the integrand 
in d2.21 i is also a polynomial so that the line integral can be exactly computed by 
means of a suitable quadrature formula. It is easy to observe that in general, due to 
the high degree of the integrand function, such quadrature formula cannot be solely 
based upon the available abscissae {c,}: one needs to introduce an additional set of 
abscissae {c\, . . . ,c,-}, distinct from the nodes {c,}, in order to make the quadrature 
formula exact: 

f 1 6 (t + Th) T VH(a(t + rh))dr = (2.8) 
Jo 

s r 

£ Pt&(to + Ci h) T VH(c>(to + ah)) + £ ftd(fo + Cih) T VH{(7{t + dih)), 

i=l z=l 

where J3,-, i = 1 , . . . , s, and Pi,i= 1 , . . . , r, denote the weights of the quadrature formula 
corresponding to the abscissae {c/} U {c;}, i.e., 




(2.9) 



Remark 3 In the case considered in the previous Remark\l\ i.e. when co = is for- 
mally added to the abscissae H2.4\l , the first product in each formula in \2.9\ ranges 
from j = to s. Moreover, also the range of {j3,} becomes i = 0, However, 
for sake of simplicity, we shall not consider this case further. 

According to ||211 . the right-hand side of ( 12. 8t is called discrete line integral, 
while the vectors 

% = a{tQ + Cih), i=l,...,r, (2.10) 

are called silent stages: they just serve to increase, as much as one likes, the degree 
of precision of the quadrature formula, but they are not to be regarded as unknowns 
since, from ( 12.6b and (12.10b . they can be expressed in terms of linear combinations of 
the fundamental stages (12.5b . 

Definition 1 The method defined by substituting the quantities in ( 12.6l l into the right- 
hand side of ( 12.8b , and by choosing the unknown coefficients {jj} in order that the 
resulting expression vanishes, is called Hamiltonian Boundary Value Method with k 
steps and degree s, in short HB VM(k,s), where k = s + r ^5]/. 



Blended HBVMs 



5 



In the sequel, we shall see that HBVMs may be expressed through different, 
though equivalent, formulations: some of them can be directly implemented in a com- 
puter program, the others being of more theoretical interest. 

Because of the equality ( 12.8b , we can apply the procedure directly to the original 
line integral appearing in the left-hand side. With this premise, by considering the 
first expansion in ( 12.61 l. the conservation property reads 

£yj [ [ Pj(z)VH(a(t + zh))dT = 0, 

.7=1 

which, as is easily checked, is certainly satisfied if we impose the following set of 
orthogonality conditions: 

Yj= f Pj(T)JVH(a(t Q + Th))dT, ./ 1 v. (2.11) 

Jo 

Then, from the second relation of (12.61 ) we obtain, by introducing the operator 

L(f;h)o(t + ch) = (2.12) 

a(t )+hV [ C Pj(x)dx f 1 Pj(z)f(a(to + zh))dx, ce [0,1], 
"Jo Jo 

that a is the eigenfunction of L(JVH;h) relative to the eigenvalue X = 1: 

a=L(TVH;h)o. (2.13) 
Definition 2 Equation ( 12.751 ) is the Master Functional Equation (MFE) defining a 

Remark 4 Some further details are in order to better elucidate the role of the MFE in 
devising our methods. First of all we observe that, by definition, the MFE intrinsically 
brings, with its polynomial solutions (J (jo + ch), the conservation property of the 
Hamiltonian function: indeed ( 12.13b is equivalent to d2.8l l under the choice d2.1 lb . 

This means that, when searching for its solutions, one should always take care 
of the precise dimension of the polynomial vector space, say V, H(y) is intended to 
belong to: the higher is V, the higher must be the number of silent stages ( and hence 
the number of steps k) to guarantee that (12.8b be satisfied. This explains the way the 
solutions of the MFE depends on k. 

It is also clear that, assuming the same kind of distribution for all the k the nodes 
(see later), ( 12.8b will be satisfied starting from a suitable number of steps k = k v on. 
This implies that, for all k > k v , HBVM(k,s) will define the same polynomial o of 
degree s, such that H(<j(tQ + h)) = H(o(to)). 

To practically compute a, we set (see ( 12.5b and ( 12.6b ) 

s 

Y i = a(t Q +c i h)=y + hY,a ij Yj, i=l,...,s, (2.14) 
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where 

aij = Pj (x)dx, i,j=l,...,s. 

Inserting ( 12.1 lb into ( 12.14b yields the final formulae which define the HBVMs class 
based upon the orthonormal basis {Pj}: 

Y i =y Q + hYa ij f 1 Pj(T)JVH(a(t + Th))dT, i = l,...,s. (2.15) 

We recall once again that we are working under the assumption (12.8b . namely 
that the Hamiltonian is a polynomial and that we are considering a sufficient num- 
ber of additional abscissae c, such that the line integral and its discrete counterpart 
do coincide. This implies that we can replace the integrals appearing in ( 12.15b by 
sums representing the associated quadrature formulae introduced in (12.81 ). without 
introducing any discretization error. 

This leads back to express the HBVM(fc,i) method in terms of the fundamental 
stages {Yj} and the silent stages {Yj} (see ( 12.101 )). By using the notation 

f(y)=JVH(y), (2.16) 

we obtain 

Y i =y + hj^a i j[j^p l Pj(c l )f(Y l ) + £fi l Pj(d l )f(Y l U , i = l,...,s. (2.17) 

7=1 V=l '=1 / 

We again stress that the silent stages % may be removed from ( 12.17b by observing 
that, for example, 

f/ = J^(f + W/, l = \,...,r, 

i=\ 

where l(t ) are the cardinal Lagrange polynomials defined on the nodes to + C(h, i = 
l,...,s. 

From the above discussion it is clear that formulae ( 12.151 ) also make sense in the 
non-polynomial case. In fact, supposing to choose the abscissae {c,} so that the sums 
in ( 12.171 ) converge to an integral as r = k — s —> °°, the resulting formula is again 
dH5]). 

Definition 3 Formula ( 12.751 ) is named °o-HBVM of degree s or HBVM(°°,s) 

This implies that HBVMs may be as well applied in the non-polynomial case 
since, in finite precision arithmetic, HBVMs are undistinguishable from their limit 
formulae ( 12.15b , when a sufficient number of silent stages is introduced. The aspect 
of having a practical exact integral, for k large enough, was already stressed in J3]|5] 
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2.1 Runge-Kutta formulation of HBVMs 

On the other hand, we emphasize that, in the non-polynomial case, (12. 1 5b becomes 
an operative method only after that a suitable strategy to approximate the integrals 
appearing in it is taken into account. In the present case, if one discretizes the Master 
Functional Equation ( 12.121 i-( 12.131 l. HBVM(£:,.s) are then obtained, essentially by ex- 
tending the discrete problem ( 12.171 ) also to the silent stages ( 12.101 ). In order to simplify 
the exposition, we shall use ( 12.161 ) and introduce the following notation: 

{r,} = { Cl }U{c,}, {co,} = {A}U{j3 ; }, 

yi = a(to+tih), fi = f(a(t Q + tih)), i=l,...,k. 

The discrete problem defining the HB VM(k, s) then becomes, 

y i = yo+hf i f'Pj(x)dxY l (»ePj(te)fe, i=l,...,k. (2.18) 

j=i Jo t=i 

By introducing the vectors 

y = {y\,...y k ) T , « = (l,...,l) r eR* 

and the matrices 

fl=dug(flh,...,flfc) ) 4^el to , (2.19) 
whose (/, 7*)th entry are given by 

(S s )ij= f'Pj{x)dx, (&s)ij=Pj(ti), (2.20) 
Jo 

we can cast the set of equations ( 12.181 ) in vector form as 

y = e ® y + h{J s & T s Q. ) g> W(y) , 

with an obvious meaning of /(y). Consequently, the method can be seen as a Runge- 
Kutta method with the following Butcher tableau: 



h 








a c 

- ' S ' 




tk 








0)1 . 





(2.21) 



Remark 5 We observe that, because of the use of an orthonormal basis, the role of 
the abscissae {c,} and of the silent abscissae {c,} is interchangeable, within the set 
{//}. This is due to the fact that all the matrices ,J? S , £P S , and Q depend on all the 
abscissae {f,}, and not on a subset of them, and they are invariant with respect to the 
choice of the fundamental abscissae {c ( -}. 

Hereafter, we shall consider a Gauss distribution of the abscissae {t\, ... ,^}, so 
that the resulting HBVM(k,s) method 0: 
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- has order 2s for all k>s; 

- is symmetric and perfectly A-stable (i.e., its stability region coincides with the 
left-half complex plane, C~ lfT4l '): 

- reduces to the Gauss-Legendre method of order 2s, when k = s; 

- exactly preserves polynomial Hamiltonian functions of degree v, provided that 




3 The Isospectral Property 



We are now going to prove a further additional result, related to the matrix appearing 
in the Butcher tableau d2.2U . corresponding to HBVM(fc,s), i.e., the matrix 

a = f.&la 6 w kxk , 



k>s, 



(3.1) 



whose rank is s. Consequently it has a (k — s)-fold zero eigenvalue. In this section, 
we are going to discuss the location of the remaining s eigenvalues of that matrix. 

Before that, we state the following preliminary result, whose proof can be found 
in El page 79]. 

Lemma 1 The eigenvalues of the matrix 



X s = 



with 



1 



(3.2) 



(3.3) 



2y/(2j +1)(2;-1) 

coincide with those of the matrix in the Butcher tableau of the Gauss-Legendre method 
of order 2s. 

We also need the following preliminary result, whose proof derives from the prop- 
erties of shifted-Legendre polynomials (see, e.g., U or the Appendix in J5])- 

Lemma 2 With reference to the matrices in K2.19\ - §2.20\ , one has 

where 



( 



-6 




\ 



with the £j defined by ( I3.3I ). 
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The following result then holds true. 

Theorem 1 (Isospectral Property of HBVMs) For all k> s and for any choice of 
the abscissae {?,■} such that the quadrature defined by the weights {©,} is exact for 
polynomials of degree 2s — 1, the nonzero eigenvalues of the matrix A in rti.il > coincide 
with those of the matrix of the Gauss-Legendre method of order 2s. 

Proof For k = s, the abscissae {f,} have to be the s Gauss-Legendre nodes, so that 
HBVM(s, s) reduces to the Gauss Legendre method of order 2s, as outlined at the end 
of Section [2] 

When k > s, from the orthonormality of the basis, see ( 12.7b . and considering that 
the quadrature with weights {©,} is exact for polynomials of degree (at least) 2s— 1, 
one easily obtains that (see d2.191 l- (l2.20b ) 

since, for all i = l,...,s, and j =l,...,s+l: 



+1, 



£ coePi(ti)Pj(te) = / Pi(t)Pj(t)dt = 8, 



By taking into account the result of Lemma|2l one then obtains: 

a^ s+1 = j? s <?fn& s+1 = J s (i s o) = & s+l x s (/,. o) = ^ s+1 (x s o) 



s+l 



-6 





V" 



-4- 



with the {£,j} defined according to ( 13.31 ). Consequently, one obtains that the columns 
of &s+l constitute a basis of an invariant (right) subspace of matrix A, so that the 
eigenvalues of (a, 0) are eigenvalues of A. In more detail, the eigenvalues of (X s 0) 
are those of X s (see ( 13.21 )) and the zero eigenvalue. Then, also in this case, the nonzero 
eigenvalues of A coincide with those of X s , i.e., with the eigenvalues of the matrix 
defining the Gauss-Legendre method of order 2s. □ 



4 Solving the discrete problem 

We shall now consider some computational aspects concerning HBVM(A:,5). In more 
details, we now show how its cost depends essentially on s, rather than on k, in the 
sense that the nonlinear system to be solved, for obtaining the discrete solution, has 
(block) dimension s. This has been already shown in J5], but here we derive the same 
result in a slightly more compact way, which will allow us to easily introduce blended 
HBVMs in the next section. 
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In order to simplify the notation, we shall fix the fundamental stages at t\ , . . . ,t s , 
since we have already seen that, due to the use of an orthonormal basis , they could be 
in principle chosen arbitrarily, among the abscissae {f,}. With this premise, we have, 
from ( 12131 

s k 

y i ^yo + hY, a ijJl 0) i p j( t ^f^ i = l,...,s. (4.1) 

7=1 f=l 

This equation is now coupled with that defining the silent stages, i.e., from ( 12.61 ) 
and d2l0b . 

S rtj 

yi=ya+Y J Yj Pj(t)dt, i = s + l,...,k. (4.2) 
y=i Jo 

Let us now partition the matrices J s , &> s G R kxs in d2.19l )-( f2~20b into J sX , 8P S \ G 
W xs and J? S 2, £? S 2 G M , containing the entries defined by the fundamental ab- 
scissae and the silent abscissae, respectively. Similarly, we partition the vector y into 
yi, containing the fundamental stages, and y2 containing the silent stages and, ac- 
cordingly, let £2\ e W xs and X2 2 G M. k ~ sxk ~ s be the diagonal matrices containing the 
corresponding entries in matrix £2. Finally, let us define the vectors y = (y[ , . . . , yj ) T , 

c=(l lfeR', and u = (1, . . . , l) r G R k ~ s . 

Consequently, we can rewrite ( 14. Il l and ( I4.21 i. as 

y, = e®yo + hAi Q ^j®h m (/(^j) < ( 4 - 3 > 

yi = u ® y + hJ s2 <g> I 2m Y, (4.4) 

respectively. The vector y can be easily retrieved from the identity ( 12.141 ). which in 
vector form reads 

yi = e ® y + hJ s \ ® 7 2m y, 

thus giving 

y 2 = (h - JfaS^e) ®y + J^J^ 1 (8 Wi 

=e wgiyo+A! ®/2 m yi, (4.5) 

in place of ( 14.41 ). where, evidently, A\ G R* - ' 5 * 1 '. By setting 

Si = ^vi ^fifli G R iXS , B 2 = ^1 ^J 2 ^2 G R sxfc_i , (4.6) 

substitution of (14.51 ) into ( 14.3b then provides, at last, the system of (block) size s to be 
actually solved: 

F(yi) = yi-e®y -h\B l ®l lm f{j{)+ ( 4 -7) 
#2 <8 /2m/ (« (Xiyo + Ai ® 7 2m yi )] = 0. 

By using the simplified Newton method for solving (14.7b . and setting 



C = B l +B 2 A l GR S 



(4.8) 
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one obtains the iteration: 



(/,®/ 2m -/zC®/ )SW = -F(y ( / 7) ) = V [ n \ (4.9) 
y (« + i) =y W + 5W) B = ,l 

where Jq is the Jacobian of f(y) evaluated at yo- Because of the result of TheoremQ] 
the following property of matrix C holds true. 

Theorem 2 The eigenvalues of matrix C in ( I4.c?| ) coincide with those of matrix ( 15.21 ), 
i.e., with the eigenvalues of the matrix of the Butcher array of the Gauss-Legendre 
method of order 2s. 

Proof Assuming, as usual for simplicity, that the fundamental stages are the first 
s ones, one has that the discrete problem 

®y + hA®I 2m f(y), 

which defines the Runge-Kutta formulation of the method, is equivalent, by virtue of 
(gl) , (1431 E2), to 

where, as usual, r = k — s0 Consequently, the eigenvalues of the matrix A defined in 
(13. U coincide with those of the matrix pencil 

Is O sxr \ ( Bi B 2 



that is 

M G a (A) n 



I s O sxr \ ( u\ _ f B\ B 2 Wu 
-Ai / r M v J lo rxJ O rxr / I v 



for some nonzero vector (u r , \ T ) T . By setting u = 0, one obtains the r zero eigenval- 
ues of the pencil. For the remaining s (nonzero) ones, it must be v = A]U, so that: 

juu= (Biu + B 2 \) = (B l u + B 2 A i u) =Cu n G a(C). □ 



Remark 6 From the result of Theorem^ it follows that the spectrum of C doesn't 
depend on the choice of the s fundamental abscissae, within the nodes {?;}. On the 
contrary, its condition number does: the latter appears to be minimized when the fun- 
damental abscissae are symmetrically distributed, and approximately evenly spaced, 
in the interval [0, 1]. As a practical "rule of thumb", the following algorithm appears 
to be almost optimal^ 

2 As observed in [19. 22], such formulation fits the framework of block BVMs. 

3 We plan to investigate this aspect further, in a forthcoming paper. 
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1. let the k abscissae {?,-} be chosen according to a Gauss-Legendre distribution of 
k nodes; 

2. then, let us consider s equidistributed nodes in (0,1), say {? i , . . . , i s }; 

3. select, as the fundamental abscissae, those nodes, among the {£,-}, which are the 
closest ones to the {tj}; 

4. define matrix C in accordingly. 

Clearly, for the above algorithm to provide a unique solution ( resulting in a symmet- 
ric choice of the fundamental abscissae), the difference k — s has to be even which, 
however, can be easily accomplished. 

In order to give evidence of the effectiveness of the above algorithm, in Figure HTTl 
we plot the condition number of matrix C = C(k, s), for s = 2, ... ,5, and k > s. As one 
can see, the condition number of C(k, s) turns out to be nicely bounded, for increasing 
values of k, which makes the implementation (that we are going to analyze in the 
next section) effective also when finite precision arithmetic is used. For comparison, 
in Figure l4~2l there is the same plot, obtained by fixing the fundamental abscissae as 
the first s ones. 



5 Blended HBVMs 

The solution of problem ( 14.91 ) is now cast into the framework of blended implicit 
methods H2ll7ll8l l9lll0lll3lll5ll23ll as below described. First of all, we observe that, 
since C is nonsingular, we can recast problem ( 14.91 ) in the equivalent form 

j{C- 1 ®h m ~hI s ®J Q ) 8 (n) = -yC- 1 ®I 2 mF(y { " ] ) = yft\ (5.1) 

where y > is a free parameter to be chosen later. Let us now introduce the weight 
( matrix) function 

G=I s (g>cp-\ <P = I 2 ,„-hYJoeR 2mx2m , (5.2) 

and the blended formulation of the system to be solved, 

MS w = [0 (/,. (8 1 2m -hC®J ) + 

(I - 9)j(C- 1 <8/ 2m -W s (g)/o)]5 W 

= 9yf i " ) + (I-e)w { 2 ) = W {n) - (5-3) 

The latter system has still the same solution as the previous ones, since it is obtained 
as the blending, with weights and (/— 9), of the two equivalent forms (14. 9t and 
(15.11 ). For iteratively solving ( 15.3b . we use the corresponding blended iteration, for- 
mally given by ElHEEHOliniEllia : 



5 (n,l+l) = 5 (n,l) _ q ( m (n,l) _ ^ ^ ^ = , 1, . 



(5.4) 



Blended HBVMs 



13 




10° ' ' 1 1 1 1 1 1 1 1 

10 20 30 40 50 60 70 80 90 100 

k 

Fig. 4.1 Condition number of the matrix C = C(k,s), for s = 2,3,4,5 and k = s,s + 1,... , 100, with the 
fundamental abscissae chosen according to the algorithm sketched in Remark[6] 




k 

Fig. 4.2 Condition number of the matrix C = C(k,s), for s = 2,3,4,5 and k = s,s + 1, . . . , 100, with the 
fundamental abscissae chosen as the first s ones. 
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Remark 7 A nonlinear variant of the iteration ( 15.41 ) can be obtained, by setting 




Y2 an d V^"'^ similarly defined, as: 

8 {n,t+l) = 8 {n,i) _ q f MS M _ ^nfi\ ; £ = 0, 1, . . . . (5.5) 

Remark 8 We emphasize that, for actually performing the iteration ( I5.2| )-( l5~?l ), as 
well as ( 15. 51 ), one has to factor only the matrix <P in ( 15.21 ), which has the same size as 
that of the continuous problem, due to the (block) diagonal structure of 9. 

We end this section by observing that the above iteration (15.4b depends on a free 
parameter y. It will be chosen in order to optimize the convergence properties of the 
iteration, according to a linear analysis of convergence, which is sketched in the next 
section. 



6 Linear analysis of convergence 

The linear analysis of convergence for the iterations d5.4l ) is carried out by considering 
the usual scalar test equation (see, e.g., ifTTl and the references therein), 

y'=Xy, 9i(A)<0. 

By setting, as usual q = hX, the two equivalent formulations (14. 9\ and ( 15. U become, 
respectively (omitting, for sake of brevity, the upper index «), 

(I s -qC)5 = Y l , r(C- 1 -qh)8 = y 2 . 

Moreover, 

e = e(q) = {\-yq)- X I Sl (6.1) 
and the blended iteration (15.41 i becomes 

5 ( ^ +1) = (/, - e(q)M(q))8 {£) + 9(q)w(q), (6.2) 

with 

M(q) = 9(q)(I s -qC) + (I s -9(q))Y{C- 1 - ql s ) , (6.3) 

¥(q) = 0(q)Wi + (Is-e{q))Wz- 

Consequently, the iteration will be convergent if and only if the spectral radius, say 
p (q), of the iteration matrix, 

Z(q)=I s -6(q)M(q), (6.4) 

is less than 1 . The set 

r = {qeC: p(q)<l} 

is the region of convergence of the iteration. The iteration is said to be (see IfTTl for 
details): 
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- A -convergent, if C CT; 

- L-convergent, if it is A-convergent and, moreover, p (q) — s- 0, as q —> °°. 
For the iteration (16.2b one verifies that (see ( 16. lb . ( 16.3b , and (16.4b ) 

Z(g) = — g^ C-^C-y/,) 2 , (6.5) 

which is the null matrix at q = and at «>. Consequently, the iteration will be A- 
convergent (and, therefore, L-convergent), provided that maximum amplification fac- 
tor, 

p* = max p(q) < 1. 

9t(,)=0 



From ( 16.5b one has thatQ 



By taking into account that 



9?( 9 )=0|(1-7^) 2 | 2 7 ' 

one then obtains that 

p* = max lf ~ , , . 
M ea(c) 2y|/i| 

For the Gauss-Legendre methods (and, then, for any matrix C having the same spec- 
trum), it can be shown that lT7l fT3l the choice 

7= iMminl = min \n\, (6.6) 

iiea{C) 

minimizes p*, which turns out to be given by 

p * = 1 - cos <p min < 1 , <p min = Arg {n min ) . (6.7) 



In Table 16.11 we list the optimal value of the parameter y, along with the corre- 
sponding maximum amplification factor p*, for various values of s, which confirm 
that the iteration ( 16.2b is L-convergent. 

Remark 9 We then conclude that the blended iteration ( 15.41 ) turns out to be L-conver- 
gent for HBVM(k, s) methods, for all s > I and k> s. 



4 Hereafter, ff(C) will denote the spectrum of matrix C. 
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Table 6.1 Optimal values )6.6t . and corresponding maximum amplification factors 46.7t . for various 
values of s. 



s 


r 


P* 


2 


0.2887 


0.1340 


3 


0.1967 


0.2765 


4 


0.1475 


0.3793 


5 


0.1173 


0.4544 


6 


0.0971 


0.5114 


7 


0.0827 


0.5561 


8 


0.0718 


0.5921 


9 


0.0635 


0.6218 


10 


0.0568 


0.6467 



7 Conclusions 

In this paper, computational aspects related to the efficient implementation of HBVM 
methods with k steps and degree s (k > s) (in short, HBVM(k,s)), have been recast 
in the framework of blended implicit methods. In more details, we have seen that 
the discrete problem generated by HBVM(k 1 s) amounts to a nonlinear system of 
(block) dimension s. Its efficient solution can be obtained by considering the blended 
formulation of the discrete problem, for which an efficient diagonal splitting can be 
easily defined. Consequently, to implement the nonlinear iteration, only one matrix 
having the same size as that of the continuous problem has to be factored. The free 
parameter, on which the blended iteration depends on, can be easily chosen, because 
of the isospectral property of HBVM(&,.s), resulting in an L-convergent iteration. 
Last, but not least, also the conditioning of the discrete problem depends only on 
s, and it appears to tend to a (nicely) bounded limit, as k grows. We plan, in the 
future, to implement HB VMs in blended formulation (in short, blended HBVMs), in 
a computational code for numerically solving Hamiltonian problems. 
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